Field-induced superconductor to insulator transition in Josephson-junction ladders 
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The superconductor to insulator transition is studied in a self-charging model for a ladder of 
Josephson-junctions in presence of an external magnetic field. Path integral Monte Carlo simulations 
of the equivalent (1 + l)-dimensional classical model are used to study the phase diagram and critical 
behavior. In addition to a superconducting (vortex-free) phase, a vortex phase can also occur for 
increasing magnetic field and small charging energy. It is found that an intervening insulating phase 
separates the superconducting from the vortex phases. Surprisingly, a finite-size scaling analysis 
shows that the field-induced superconducting to insulator transition is in the KT universality class 
even tough the external field breaks time-reversal symmetry. 
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Superconductor to insulator transitions in Joseph- 
son junction arrays have attracted considerable 
interes1pi*2i2iii£. Such arrays can currently be fabri- 
cated in any desired geometry both in one and two 
dimensions 4 - with well-controlled parameters. When 
charging effects due to the small capacitance of the grains 
or junctions dominate, strong quantum fluctuations of 
the phase of the superconducting order parameter 
may drive the system into an insulating phase at zero 
temperature leading to a superconductor to insulator 
transition as a function of charging energy or external 
magnetic field. In two dimensions, the universality class 
of these transitions have already been investigated in 
detail numerically 6 , both in relation to experiments^ 
and theoretical predictions^^. Nevertheless, there are 
also remarkable quantum phase transition scenarios^ 
that can take place in a Josephson-junction ladder in 
a magnetic field as in Fig. 1 which have not been 
investigated experimentally or numerically. Inter- 
estingly, a Josephson-junction ladder provides the 
simplest one-dimensional version of the array allowing 
the study of commensurability effects due to the flux 
latticetlLUJJdy^i or charge frustration 14 , in presence of 
quantum fluctuations. 

An important aspect of the Josephson-junction lad- 
der in Fig. 1 is the approximate relation to the quan- 
tum sine-Gordon modeli^, as first shown by Kardar— . 
For increasing magnetic field there is a transition from 
a (vortex free) superconducting phase to a vortex 
phase where flux penetrates the ladder, in absence of 
quantum fluctuations. This transition is the analog 
of the commensurate-incommensurate transition 1 ^ de- 
scribed by the sine-Gordon model. For small fields, the 
phases in different branches of the ladder are locked to 
each other while in the vortex state exponentially in- 
teracting kinks (vortices) appear that unlock the phases 
leading to a one-dimensional vortex lattice. Inclusion of 
quantum fluctuations due to charging effects leads to the 
interesting prediction that an insulating phase should oc- 
cur in the vicinity to this transition and therefore a direct 
transition from the (vortex free) superconducting phase 



to the vortex lattice phase is not possible even for small 
charging energies. So far, this remarkable collective ef- 
fect resulting from small capacitances in the ladder has 
not been observed experimentally or even numerically. 
The quantum critical behavior is also of great interest. 
While at zero field the superconductor to insulator tran- 
sition is in the Kosterlitz-Thouless universality classii*i£, 
the field-induced superconductor to insulator transition 
could be in principle in a different universality class since 
the magnetic field breaks the time reversal symmetry 
and the commensurate-incommensurate transition in the 
sine-Gordon model is in a distinct universality classi£*i&. 

In this work, we use path integral Monte Carlo (MC) 
simulations of the equivalent (I + I )-dimcnsional classi- 
cal model to study the phase diagram and critical behav- 
ior of the Josephson-junction ladder. The exchange MC 
method (parallel tempering) 1 ^ is used to determine more 
accurately the scaling behavior. It is found that an in- 
tervening insulating phase separates the superconducting 
from the vortex phases, as shown in Fig. 2, in agreement 
with Kardar's prediction 1 ^ Thus, for increasing field a 
single transition to an insulator occurs for charging en- 
ergies in the range 9b < 9 < 9a while for lower charging 
energies the superconductor to insulator transition is fol- 
lowed by a transition into a vortex phase. Surprisingly, 
a finite-size scaling analysis shows that the field-induced 
superconducting to insulator transition AC is in the KT 
universality class. 

We consider a periodic Josephson-junction ladder as 
indicated in Fig. 1, where charging effects are domi- 
nated by the capacitance to ground of each grainii, and 
described by the Hamiltoniar 
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The first term in Eq. Q describes quantum fluctuations 
induced by the charging energy E c = 4e 2 /C of a non- 
neutral superconducting grain located at site r, where 
e is the electronic charge and C is the effective capaci- 
tance to the ground of the grain, while the second term is 
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FIG. 1: Schematic representation of a periodic Josephson- 
junction ladder. Superconducting grains (squares) with 
charging energy E c are coupled to the nearest neighbors by 
intra (E x ) and inter-chain (E y ) Josephson couplings (crosses). 
The frustration / corresponds to the magnetic flux through 
the elementary cell, / = $/$o, in units of the flux quantum 



the usual Josephson-junction coupling between nearest- 
neighbor grains described by phase variables 9 r . The 
effect of the applied magnetic field appears through the 

bond variables A rr i = (2n/& ) A ■ dr, where A is the 
vector potential due to the external magnetic field B and 
the gauge-invariant sum around an elementary cell of the 
ladder is given by Y^rr' = w ' tn / = ^/^o defining 
the frustration parameter. It is sufficient to consider / 
in the range [0, 0.5]. 

In order to study the effects of the magnetic field on 
the superconductor to insulator transition and the crit- 
ical behavior, it is convenient to use an imaginary-time 
path-integral formulation of the model£. In this formula- 
tion, the ID quantum problem of Eq. (1) maps into a 2D 
classical statistical mechanics problem with the extra di- 
mension corresponding to the imaginary-time direction. 
The time axis r can be discretized in slices At and the 
ground state energy of the quantum model corresponds 
to the reduced free energy of the classical model, per unit 
length in the imaginary time direction. After re-scaling 
the time slices appropriately in order to get space-time 
isotropic couplings, the resulting classical partition func- 
tion is given by Z = Tr^gye~ H where the reduced classi- 
cal Hamiltonian can be defined as 
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In the above equation, 9 and 9' denote the phases in the 
upper and lower branches in Fig. (1), j and r label the 
spatial and time directions, respectively, and the ratio 
of the charging energy to the Josephson coupling g — 
(E c / Ex) 1 / 2 plays the role of an effective temperature in 
the 2D classical model. The equation © corresponds to 
a gauge choice where the vector potential A is parallel 
to ladder, with opposite values in the upper and lower 
branches. 

We carry out MC simulations using the 2D classical 
Hamiltonian in Eq. [21 regarding g as a temperature- like 
parameter and employing the standard Metropolis's al- 
gorithm to generate the equilibrium distribution. Near 




FIG. 2: Phase diagram of the ladder obtained from path inte- 
gral MC simulations showing the (vortex-free) superconduct- 
ing phase (S), the insulating phase (I) and the vortex phase 
(V), as a function of frustration / and ratio of charging to 
Josephson energy g = (E c /E x ) 1 '' 2 for E y /E x = 2. The criti- 
cal field f c at g — was obtained separately, by considering 
the ground state of the ladder in absence of quantum fluctu- 
ations (E c = 0). The lines are just guides to the eye. 



the superconductor to insulator transition, we use the ex- 
change MC method (parallel tempering}— to determine 
more accurately the scaling behavior of various physical 
quantities. In fact, this method is known to reduce sig- 
nificantly the critical slowing down near the transition. 
In this method, many replicas of the system with differ- 
ent couplings j in a range above and below the critical 
point are simulated simultaneously and the correspond- 
ing configurations are allowed to be exchanged with a 
probability distribution satisfying detailed balance. Sim- 
ulations are performed in system sizes with equal spatial 
and time linear length L. This choice of the aspect ratio 
of the system assumes implicitly that the dynamic critical 
exponent z characterizing the superconductor to insula- 
tor transition is close to z ~ 1. In general, a quantum 
phase transition is characterized by intrinsic anisotropic 
scaling with different diverging correlation lengths £ and 
£ r in the spatial and time directions^, respectively, re- 
lated by the dynamic exponent z as £ r oc £ z . Our choice 
is justified by the analysis discussed below, showing that 
the scaling behavior is in fact consistent with z = 1. 

To locate the superconductor to insulator transition 
as a function of charging energy and magnetic field, we 
have performed MC calculations of the helicity modulus 
(phase stiffness), p and p T , in the spatial and time di- 
rections. In the superconducting phase these quantities 
should be finite, reflecting the existence of phase coher- 
ence, while in the insulating phase they should vanish in 
the thermodynamic limit. Fig. 3a shows the behavior 
of the helicity modulus as a function of frustration for a 
fixed value of g, below the critical value corresponding to 
/ = 1/2 frustration in Fig. 2, g B ~ 0.96. The helicity 
modulus remains finite at small fields and at fields close 
to / = 1/2 where a superconducting phase is expected-^ 
but at intermediate fields reach small values which de- 
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FIG. 3: Phase stiffness in spatial p and imaginary time p T di- 
rections (a) and energy fluctuation C — (< H 2 > — < H > 2 
)/L 2 of the equivalent classical model (Eq. @), as a function 
of frustration /, for a constant value of the ratio g — 0.9. 
The double peaks in C indicate two transitions for increas- 
ing / with an intervening insulating phase corresponding to 
vanishing phase stiffness. 



creases with increasing system sizes. Fig. 3b also show 
the behavior of the energy fluctuation (specific heat of 
the classical model) C as a function of frustration show- 
ing two peaks near the two critical fields defining the 
region where the helicity modulus reach small values. 
The behavior in Figs 3a and 3b indicate that there are 
two superconducting to insulator transitions for increas- 
ing / with an intervening insulating phase corresponding 
to vanishing phase stiffness for large system sizes L. Per- 
forming the same calculations for different couplings g, 
we have constructed the phase diagram in Fig. 2 for 
E y /E x — 2 showing three different phases separated by 
the transition lines AC and BC with the insulating phase 
extending between the superconducting vortex-free phase 
S and vortex phase V for g < gs- The critical field f c 
at g — was obtained separately, by considering the 
ground state of the ladder in absence of quantum fluctu- 
ations [E c = 0). Using simulated annealing to obtain the 
ground state, f c was determined as the value of / where 
flux first penetrates the ladder for increasing /. Note 
that, due to the long runs required to reach proper equi- 
libration, the lowest nonzero value of g accessible within 
the present MC calculations is g ~ \gB- Therefore, this 
result does not rule out the possibility that the two tran- 
sitions lines in Fig. 2 merge at a lower value of g into a 



single line ending at f c in which case a direct transition 
from S to V phases is possible for sufficiently small g. We 
have performed additional calculations with various ra- 
tios E y /E x and found that the insulating phase is more 
clearly visible for E y /E x > 1.5 while at lower values, the 
transition lines AC and BC are difficult to resolve at low 
values of g. 

The intervening insulating phase between the super- 
conducting phase at low fields and vortex phase at higher 
fields was predicted analytically in an earlier work by 
Kardar 10 using an approximate relation of the Hamilto- 
nian of Eq. with Gaussian and sine-Gordon mod- 
els describing fluctuations of the phase variables cj) = 
(9' + 6)/2 and ip = (9' - 9)/2, respectively, where 9' 
and 9 denote the phases in the upper and lower branches 
of the ladder in Fig. (1). This approximation allows 
one to calculate the power-law correlation functions for 9 
and 9' , C(r) oc 1/r'', describing the superconducting and 
vortex phases, using the exact results for the sine-Gordon 
model 15 . However, the assumed Gaussian approximation 
leaves out space-time vortices which are responsible for 
the insulating phase where the correlation function de- 
cays exponentially. The insulating phase was then stud- 
ied by adding the effects of space-time vortices "by hand" 
through the usual KT criterion r\ < 1/4 to determine 
the bound of stability of the power-law correlated phase. 
The results showed that the insulating phase expected 
for large g, extends all the way to g = at f c . The phase 
diagram in Fig. 2 is consistent with this prediction but 
as mentioned above a direct transition between S and V 
phases at sufficiently small g can not be ruled out within 
the present calculations. 

Although a KT-like stability criterion was used in Kar- 
dar's workiS to determine approximately the location of 
the transitions to the insulating phase in the phase di- 
agram, the nature of the transition can not be deter- 
mined by such argument because it is based in decou- 
pled Gaussian and sine-Gordon models. Like in a similar 
analysis employing the KT criterion for the sine-Gordon 
modeliS, applied to adsorbed layers, which also predict 
an intervening disordered phase near the commensurate- 
incommensurate transition for low-order commensurate 
phases, the universality class of the commensurate to 
disordered phase has to be determined by other methods 
based on symmetry considerations, mapping to a solvable 
model or numerical simulations. 

The universality class of the field-induced supercon- 
ductor to insulator transition along the line AC in Fig. 
2 is of particular interest. At zero magnetic field (/ = 0) 
the transition as a function of charging energy is in 
the KT universality classiiii while at nonzero magnetic 
fields in principle a different universality class is possible 
since the magnetic field breaks the time-reversal sym- 
metry. In addition, the relation of the ladder model 
of Eq. (1) with the quantum sine-Gordon chainiS also 
suggests the possibility of a universality class different 
from the KT transition at zero field since the super- 
conductor to insulator transition in the ladder is driven 
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FIG. 4: Behavior of the scaled correlation length £ T /L (a) 
and £/L (b) across the superconductor to insulator transition 
line AC in Fig. 2 at a fixed nonzero frustration / = 0.2, for 
different systems sizes L. 
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FIG. 5: Behavior of the phase stiffness p T (a) and p (b) across 
the superconductor to insulator transition line AC in Fig. 2 
at a fixed nonzero frustration / = 0.2, for different systems 
sizes L. 



by the commensurate-incommensurate transition in the 
sine-Gordon model which is in a distinct universality 
classi£*i&. In order to study the critical behavior of 
the transition line AC in detail we have investigated the 
finite-size behavior of the phase correlation lengths £ and 
£ r in addition to the helicity modulus p and p T near 
the transition line AC. The correlation lengths £(L) and 
£ T (L) in the finite-size system can be obtained from a 
second moment calculation using the correlation function 
asSfi 

eW = 9 . * [£(o)Ag(fc)-i] 1/2 (3) 

where S(k) is given by the Fourier transform 

S(k) =Y,< e l ^- e ^+* > e lkx (4) 

X 

The wavector k — 2n/L is the smallest value allowed in 
the finite system. Similar expressions are used for the 
correlation length £,{L) T in the time direction. 

Fig. 4 shows the finite-size behavior of the correlation 
length scaled by the system size L near the transition 
line AC in Fig. 2, as function of g, for a fixed value 



of frustration / = 0.2 < fc- In the insulating phase 
for g > g c , the correlation length is finite and there- 
fore £/L and £ T /L decrease with system size and should 
vanish in the limit L — > oo while at the transition and 
in the superconducting phase, where the system is crit- 
ical, £/£ should approach a finite value depending of g. 
The behavior in Fig. 4a and 4b, showing that both £/L 
and £ T /L merge for low values of g at approximately the 
same critical value g c ~ 1.3 is consistent with a transi- 
tion in the KT universality class where the correlation 
length diverges exponentially £ oc e b /' 9 ~ 9c . In gen- 
eral, for a power-law correlation length £ oc \g — <? c | -1 ^ 2 , 
the dimcnsionlcss quantity £/L is expected to cross at 
the transition for different system sizes2& and satisfy the 
scaling relation £/L = G((g — g^L 1 /") near the transi- 
tion. Then, the slopes s(L) of the curves £/L for dif- 
ferent L are determined by the critical exponent v as 
s(L) — oc L 1 /", evaluated near g c . The ex- 

ponential correlation length of the KT transition corre- 
sponds to v — > oo which implies that s(L) should be 
L-independent and therefore the curves £/L for different 
system sizes should merge at the critical point, as in fact 
observed in Fig. 4. 

The finite-size behavior of the helicity modulus p and 
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FIG. 6: Behavior of the quadratic average of the phase stiff- 
ness (pp T ) 1//2 across the superconductor to insulator transi- 
tion line AC in Fig. 2 at a fixed nonzero frustration / = 0.2, 
for different systems sizes L. Dotted line corresponds to the 
universal jump prediction for a KT transition. 



p T is also consistent with a KT transition. As shown in 
Fig. 5, p and p T also merge for decreasing g at approxi- 
mately the same critical g c . At the transition and in the 
superconducting phase, where the system is critical, these 
quantities should scale as2*& p oc L 2 ~ d ~ z and p T oc L z ~ d , 
in d spatial dimensions. Since p T and p in Fig. 5a and 5b 
approach finite values independent of L in this regime, 
the dynamic exponent is z = 1 for the present case where 
d = 1. This value of the dynamic exponent z means that 
the line of Gaussian fixed points can be made isotropic 
even though the phase stiffness, p and p T , are different 
in the spatial and time directions. In fact, by a suitable 
rescaling of the coordinates in the spatial and time direc- 
tions, the stiffness of the renormalized isotropic Gaussian 
model can be written as p = (p pr) 1 ! 2 . If the transition 
is in the KT universality class then a universal jumpSi is 
expected for p as a function of g at the transition, given 
by p — ~g. Fig. 6 shows that the finite-size behavior of p 
appears consistent with the expected universal jump for 
the phase stiffness. For each system size, the intersection 
of the data with the dotted line representing the univer- 
sal jump prediction provides an upper bound estimate 
of the critical coupling. A more quantitative determina- 
tion of the critical coupling and the corresponding jump 
in the phase stiffness can be obtained by the Weber and 
Minnhagen's scaling analysis, which has been introduced 
before as an accurate method of locating the critical cou- 
pling and the universal jump of the KT transition for the 
ordinary XY models. The phase stiffness at the transi- 
tion should have a logarithmic correction given by22. 
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FIG. 7: Numerical fit of the phase stiffness p according to We- 
ber and Minnhagen's scaling relation with two free parameters 
for system sizes L = 6,8, 12, 16, 18, 24, at different couplings g 
across across the superconductor to insulator transition, (a) 
The estimate of the jump pn/g c from the fitting parameter 
and (b) error estimate of the fit. The location of the minimum 
in (b) determines the critical coupling g c . Dotted line in (a) 
is the value of the universal jump for a KT transition. 
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where C is an undetermined constant and the universal 
jump implies pu = —g c . Using data for different sys- 
tem sizes, the best numerical fit according to this scaling 
form performed for different values of the coupling g can 
then be used to locate the transition and determine the 
jump. Treating both pn and C as free parameters, we 
have performed this fit using data of the quadratic av- 
eraged phase stiffness g for smaller system sizes L < 24, 
where the data is more accurate. As shown in Fig. 7, the 
value of the jump near the location of the minimum in the 
fitting error for different couplings g, is indeed consistent 
with the value expected for a KT universality class. 

In conclusion, we have studied the superconductor to 
insulator transition in a ladder of Josephson-junctions 
under an applied magnetic field by path integral MC 
methods. We found clear numerical evidence in the phase 
diagram of an intervening insulating phase between the 
superconducting and vortex phases, in good agreement 
with an earlier prediction by Kardar— . In addition, we 
find that the field-induced superconducting to insulator 
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transition is in the KT universality class. In the vor- 
tex phase, commensurability of the vortex lattice and 
the ladde»ii*i^ will strongly depend on the frustration 
parameter /. Since Josephson-junction arrays can cur- 
rently be fabricated in any desired geometry and with 
well-controlled parameters, our numerical evidence of the 
intermediate insulating phase and the universality class 



of the transition should be interesting enough to motivate 
experiments in these intriguing systems. 
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